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Abstract 


The  Defense  Nuclear  Agency  (DNA)  is  improving  the  military’s  capability  to 
forecast  dosage  and  hazard  levels  due  to  release  of  chemical,  biological  and  nuclear 
agents.  During  Operation  DESERT  STORM  the  military  realized  the  need  for  models  to 
predict  risk  levels  for  military  personnel  assigned  proximate  to  missile  attacks.  One 
project  associated  with  this  is  the  continuing  development  of  the  Operational  Multiscale 
Environment  Model  with  Grid  Adaptivity  (OMEGA).  DNA  has  sponsored  AFIT  to 
validate  OMEGA  with  focus  on  incorporating  weather  data  obtained  from  Air  Force 
Combat  Climatology  Center  for  the  period  of  the  Chernobyl  Nuclear  Accident.  The 
physies  of  the  model  is  tested  using  National  Weather  Service  Medium  Range  Forecast 
data  by  comparing  predicted  wind  fields  for  three  weather  stations  with  analysis  maps. 
The  model  is  further  tested  using  the  data  generated  at  Air  Force  Global  Weather  Central 
for  the  first  three  days  following  the  release  at  the  Chernobyl  Nuclear  Plant.  A  user- 
defined  source  term  was  developed  to  simulate  the  release  of  radionuclides  from  the 
plant.  Analysis  from  paired  t-tests  shows  statistically  how  well  OMEGA  predicts  wind 
fields.  The  results  show  qualitatively  the  promise  of  OMEGA  to  meet  the  needs  of  the 
Defense  Nuclear  Agency  as  the  model  is  still  under  development. 
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Validation  of  Operational  Multiscale  Environment  Model 
With  Grid  Adaptivity  (OMEGA) 


1.  Introduction 

The  Defense  Nuclear  Agency  is  interested  in  improving  its  capability  to  forecast 
dosage  and  hazard  levels  due  to  release  of  chemical,  biological,  or  nuclear  agents. 
Accidents  of  the  magnitude  of  the  Chernobyl  Nuclear  Accident,  explosions  like  Mount 
Saint  Helens,  and  the  chemical  spill  in  India  show  the  need  for  forecast  models  to  predict 
the  transport  of  agents.  During  Operation  DESERT  STORM  Scud  missile  attacks  by  Iraq 
provided  such  a  challenge  to  the  military.  The  ability  to  predict  transport  of  these  agents 
after  release  is  vital  to  the  survival  of  our  military  and  therefore  mission  enhancing. 
During  a  crisis  a  model  could  be  running  continually  with  real  or  simulated  sources.  This 
would  allow  full  modeling  of  conditions  within  the  domain  of  interest  before  sources  are 
introduced  and  produce  better  results. 

Traditionally,  atmospheric/mass  transport  models  use  rectangular  grids  or  polar 
stereographs  for  terrain  with  coarse  grid  resolution.  The  OMEGA  model  grid  is 
unstructured  horizontally;  it  currently  adapts  to  underlying  surface  features.  The  model 
has  a  variable  horizontal  grid  resolution  that  can  range  from  1  km  to  100  km  and  a 
vertical  resolution  ranging  from  a  few  meters  to  1  km. 

Advances  in  models  resolve  local  perturbations  on  the  larger  scale  wind  field. 
These  techniques  require  considering  physical  variables  and  processes  affecting  the  flow. 
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such  as  topography,  land  use,  land/water  interfaces,  vegetation,  soil  moisture,  surface 
moisture,  and  energy  budgets. 

OMEGA  uses  a  non-hydrostatic  equation  set  to  describe  the  dynamics. 

Turbulence  and  its  effects  are  modeled  by  a  first-order  boundary  layer  model  or  by  a 
second  order  turbulent  kinetic  energy  model  that  uses  a  prognostic  equation  for  the 
turbulent  kinetic  energy.  Cloud  formation,  growth,  and  precipitation  processes  are 
simulated  using  bulk-water  parameterization  schemes.  OMEGA  uses  a  modified  Kuo 
scheme  to  calculate  the  vertical  redistribution  of  heat  and  water  vapor  in  columns  where 
the  potential  for  deep  convection  exists.  A  convective  parameterization  scheme  is  used  in 
regions  of  coarse  resolution  where  the  potential  exists  for  vertical  accelerations  to  occur. 
The  model  also  includes  a  radiation  package  to  approximate  the  effects  of  atmosphere  and 
clouds  on  radiation. 

The  model  incorporates  an  integral  Aerosol  Diffusion  Model  into  the  design.  This 
Lagrangian  module  can  follow  massive  and  massless  aerosol,  with  user  defined  source 
and  aerosol  characterization.  Aerosols  can  be  treated  as  discrete  particles  or  as  centroids 
of  puffs  whose  dispersion  is  dependent  on  the  ambient  conditions.  OMEGA  also  has  the 
capability  to  transport  Eulerian  tracers  (massless)  to  simulate  continuous  sources  of 
gaseous  material. 

For  more  information  about  OMEGA  the  interested  reader  may  refer  to  the 
Worldwide  Web  for  on-line  documentation  describing  the  module  and  defining  the 
variables  used  for  OMEGA  at  the  following  address;  http;//tornado.saic.com. 
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Problem  Statement: 


An  alternative  to  Gaussian  plume  models— greatly  limited  with  their  simplified 
assumptions— is  needed  to  predict  atmospheric  dispersion.  While  OMEGA  appears  to 
have  many  advantages  over  current  models  and  the  Gaussian  models  still  used  by  the  U.S. 
Environmental  Protection  Agency,  it  requires  validation.  OMEGA  must  show  the  ability 
to  accurately  model  wind  flow  and  dispersion  in  various  weather,  terrain,  and  release 
conditions.  It  is  recognized  that  no  single  scenario  can  be  used  to  validate  a  model; 
limitations  of  the  model  will  be  investigated  as  well  as  strengths.  Other  runs  would  be 
necessitated  under  different  conditions  and  different  data  sets. 

The  goals  of  the  research  are  twofold.  The  first  goal  is  to  determine  how  well 
OMEGA  can  simulate  the  Chernobyl  Nuclear  Accident.  The  ability  of  the  model  to 
ingest  High  Resolution  Analysis  System  data  and  predict  deposition  of  Cs-137  will  be 
tested;  also,  other  conditions  of  the  model  such  as  length  of  simulation  and  treatment  of 
wind  fields  will  be  analyzed.  The  second  goal  will  be  to  describe  OMEGA’S 
performance  with  sample  runs  over  the  United  States  using  Medium  Range  Forecast  data 
provided  by  the  National  Weather  Service  by  comparing  predicted  wind  fields  with 
observed  weather  data.  Accomplishment  of  these  goals  will  partially  assess  the  validity 
of  OMEGA  to  meet  the  needs  of  the  Defense  Nuclear  Agency. 
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2.  Literature  Review 


2.1  Introduction 

The  demand  for  atmospheric  dispersion  models  comes  from  the  need  for  both 
diagnostic  and  prognostic  analysis  concerning  the  space  and  time  evolution  of  pollutant 
dispersal  as  a  basis  for  public  authorities  to  act  concerning  human  health  and  the 
environment.  Dispersion  models  play  an  important  role,  especially  during  the  first  stage 
of  a  pollution  episode  when  field  measurements  are  missing  or  insufficient  (Desiato, 
1992:2805-2806).  For  example,  the  U.S.  Environmental  Protection  Agency  (EPA)  needs 
air  quality  simulation  models  that  will  estimate  short-term  (hours  to  a  day)  to  long-term  (a 
few  days  to  a  year)  impact  of  source  configurations. 

The  Chernobyl  Nuclear  Accident  well  portrayed  the  need  for  such  models.  On 
26  April  1986,  31  people  were  killed  and  debris  was  scattered  for  miles  as  a  cloud  with 
radioactive  particles  and  gases  was  released  from  the  Chernobyl  Nuclear  Power  Station 
(Jagger,  1991:1-10).  The  cloud  dispersed  and  radioactive  material  was  deposited 
throughout  the  Northern  Hemisphere.  By  means  of  integration  of  environmental  data,  it 
is  estimated  that  10^’  becquerels  of  Cs-137  were  released  during  and  subsequent  to  the 
accident  (Anspaugh  et  al.,  1988: 1513).  This  corresponds  to  less  than  2  kg  of  Cs-137.  It 
is  estimated  17,000  will  die  prematurely  over  the  next  50  years  because  of  radiation, 
representing  increased  deaths  of  about  0.01%.  Dosages  in  the  Chernobyl  area  averaged 
50  rems.  Of  Hiroshima  survivors,  about  4%  received  dosages  greater  than  100  rems 
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(Jagger,  1991:124-127,  222-223).  For  comparison,  the  amount  of  Cs-137  activity 
released  from  weapon  tests  of  the  late  1950s  and  1960s  is  1.5  x  lO'*  Bq  and  from  the 
1957  Windscale  accident  4.4  x  10‘^  Bq  (Smith  and  Clark,  1989:4). 

The  initial  information  on  the  Chernobyl  release  and  its  time  variation  was  made 
available  by  the  Soviet  authorities  at  the  International  Atomic  Energy  Agency  (IAEA) 
experts  meeting  on  25-29  August  1986  in  Vienna.  This  included  estimates  of  the  total 
quantities  released  of  a  range  of  individual  nuclides,  corrected  for  decay  until  6  May  1986 
when  the  reactor  was  finally  sealed,  and  the  magnitude  of  their  release  during  the  first 
day.  At  the  1st  Steering  Committee  meeting  of  the  Atmospheric  Transport  Model 
Evaluation  Study  (ATMES)  exercise.  Dr.  Petrov  (USSR)  provided  additional  information 
giving  the  daily  release  of  Cs-137  with  an  estimate  of  the  effective  release  height  during 
the  time  after  the  initial  release  based  on  radiological  information  inside  and  outside  the 
30  km  zone  around  the  reactor,  and  model  calculations  of  the  material  deposited 
(Klug,  1992:2).  The  release  of  radioactivity  from  the  Chernobyl  reactor  occurred  as  a 
continuous  process  over  a  nine-day  period.  The  initial  release  immediately  accounted 
only  for  25%  of  the  total  radionuclides  injected  into  the  atmosphere.  After  extensive 
analysis  of  the  observations  of  radioactivity  in  the  surroundings  of  the  Chernobyl  reactor, 
it  was  concluded  that  there  were  four  distinct  stages  in  which  radionuclides  were  released. 
The  initial  stage  was  a  very  intensive  release.  During  the  following  five  days,  the  release 
rates  declined  steadily  to  a  minimum  value  six  times  lower  than  the  initial  release  rate. 
This  decline  in  the  amount  of  material  injected  into  the  atmosphere  was  then  followed  by 
a  four-day  period  of  increase  to  reach  a  level  of  70%  of  the  initial  release  rate.  The  last 
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stage  of  the  process  took  place  9  days  after  the  accident  when  the  release  was  reduced  to 
1%  of  its  initial  value. 

The  amount  released  was  estimated  on  the  basis  of  radiation  measurements  and 
various  technical  analyses  of  samples  taken  from  the  vicinity  of  the  Chernobyl  reactor.  It 
was  concluded  that  about  10‘^  to  2  x  lO'*  Bq  were  discharged  to  the  atmosphere  from  the 
core  of  the  Chernobyl  reactor.  The  observation  indicated  that  about  20%  of  volatile 
radionuclides  like  Iodine,  Cesium,  and  Tellurium  were  expelled  from  the  damaged 
reactor.  It  is  commonly  accepted  that  the  error  of  estimate  of  the  Chernobyl  release  is 
+50%.  While  experimental  data  indicate  intermittent  changes  in  the  composition  of  the 
release,  as  a  first  approximation  it  is  assumed  that  each  isotope  varied  in  time 
proportionally  to  the  total  release  (Pudykiewicz,  1990:213-225). 

The  atmospheric  release  of  radionuclides  following  the  Chernobyl  reactor 
accident  gave  rise  to  two  distinct  but  related  activities  in  the  field  of  atmospheric 
dispersion  modeling:  from  a  purely  scientific  view,  it  provided  an  opportunity  to  test 
long-range  models  developed  for  various  purposes  using  radiological  data  collected 
through  several  days  after  the  accident  over  Europe;  it  also  outlined  the  need  for 
operational  models  capable  of  describing  in  real  time  the  long  distance  transport  of  toxic 
materials  accidentally  released  into  the  atmosphere. 

The  need  for  a  computational  system  to  predict  environmental  consequences  of  a 
large-scale  nuclear  accident  has  been  recognized  by  many  countries,  especially  in  Europe, 
after  the  Chernobyl  accident.  Many  international  agencies  such  as  IAEA  and  the  World 
Meteorological  Organization  (WMO)  recognized  the  need  also  (Ishikawa,  1994:969). 
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Following  the  Chernobyl  accident  the  International  Nuclear  Safety  Advisory  Group 
recommended  IAEA  should,  in  collaboration  with  WMO,  review  and  intercalibrate 
models  of  atmospheric  transport  of  radionuclides  over  short  and  long  distances  for 
radionuclide  deposition  on  terrestrial  surfaces  and  establish  a  data  base  for  validation 
studies  on  models.  Measurements  of  radioactivity  performed  following  the  Chernobyl 
accident  provide  a  unique  database  for  extensive  verification  studies  of  atmospheric 
models.  These  international  agencies  sponsored  ATMES.  The  interested  reader  may 
refer  to  the  ATMES  exercise  where  22  models  were  compared  and  ranked  according  to 
performance  in  several  categories  (Klug,  1992:1-10). 

2.2  Elements  of  Atmospheric  Transport  Models 

General  requirements  of  any  atmospheric  dispersion  model  that  calculates  real¬ 
time  estimates  of  pollutant  air  concentration  and  deposition  fields  are:  it  must  use 
meteorological  input  data  available  in  real  time;  it  must  be  flexible  in  the  range  of 
sources,  pollutants,  and  domains  that  it  can  consider;  and  it  must  produce  results  easily 
interpreted  (Hummel  et  al.,  1990:421-426). 

One  prominent  example  is  a  Gaussian  emission  model,  probably  the  most  widely 
used  for  estimating  pollutant  dispersion  (Carrascal  et  al.,  1993: 147).  An  emission  model 
can  solve  for  plume  rise  to  estimate  effective  height  of  plume  from  a  single  source.  It  can 
estimate  concentrations  for  receptors.  Source,  height,  emission  type,  and  time  are  needed. 
The  meteorology  data  needed  includes  wind  characteristics,  Pasquill  Stability  Classes, 
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temperature,  mixing  height,  atmospheric  modeling  of  the  chemistry  and  physics,  etc. 
Receptor  coordinates  and  height  are  needed.  Then,  an  estimate  of  concentration  can  be 
made  (Turner,  1994: 1-8).  The  Gaussian  Plume  Model  will  be  discussed  in  a  later  section. 

Wind  data  should  be  frequent  enough  so  that  the  distance  of  air  mass  advection 
for  one  step  is  less  than  the  distance  between  grid  points.  This  requirement  based  on  the 
simple  relation  between  spatial  variation  of  the  wind  and  the  length  of  advection  of  one 
time  step.  Advection  beyond  the  next  grid  point  will  not  follow  the  wind  field  if  the  wind 
field  changes  in  the  distance  between  grid  points  (unless  the  spatial  variation  of  the  wind 

is  not  considered  in  the  departure  point  calculation). 

Complex  terrain  models  are  needed  for  point  sources  to  accurately  calculate  long- 
range  transport,  topographic  effects,  and  aerodynamic  downwash.  Estimation  of  pollutant 
concentrations  in  areas  of  complex  terrain  is  especially  difficult.  Under  unstable 
conditions  elevated  pollutant  plumes  have  the  tendency  to  rise  over  terrain  obstructions, 
although  they  may  pass  near  the  crest.  Under  stable  conditions,  plumes  tend  to  alter  their 
paths  to  flow  around  obstructions. 

Land-water  interface  problems  require  consideration  of  changes  in  roughness  and 
the  vertical  thermal  structure.  For  instance,  during  spring  and  summer  when  bodies  of 
water  tend  to  be  colder  than  land  surfaces  in  the  daytime,  the  potential  exists  for  elevated 
sources  in  shoreline  zones  to  fumigate  for  longer  than  the  typical  half-hour  periods  of  the 
break-up  of  nocturnal  radiation  inversions  (Turner,  1979:502-519). 

In  the  simulation  of  regional  or  hemispheric  transport,  the  curvature  of  the  earth 
should  be  included.  A  straightforward  way  is  to  describe  governing  equations  in 


8 


spherical  coordinates.  An  alternative  way  is  to  project  the  ‘spherical’  computational 
domain  onto  ‘flat’  computational  space  using  a  conformal  map  projection  (Ishikawa, 
1994:969-978).  OMEGA  structures  the  grid  such  that  the  vertices  of  the  grid  triangles  lie 
on  the  same  radius  of  the  Earth.  Each  grid  volume  is  a  truncated  triangular  pyramid  with 
the  apex  of  the  pyramid  at  the  center  of  the  Earth.  Thus,  it  uses  a  conformal  map 
projection. 

Another  factor  is  the  type  of  source.  Quantities  of  pollutants  are  released  from 
natural  sources,  industrial  sources,  mobile  sources,  and  in  many  combinations.  Some 
require  models  accounting  for  atmospheric  reactions  and  transformations  of  pollutants 
other  than  in  very  simplistic  ways. 

Models  predicting  concentration,  whether  air  or  ground  level,  must  include 
aerosol  loss  mechanisms  such  as  wet  and  dry  deposition  and  settling.  These  mechanisms 
are  often  referred  to  as  scavenging.  Dry  deposition  occurs  as  a  result  of  gravitational 
settling.  Rainout  occurs  when  particles  are  incorporated  into  droplets  while  aerosols  are 
in  clouds. 

2.3  Complications  of  Modelling 

Modelling  atmospheric  dispersion  with  the  necessary  elements  described  in 
previous  sections  is  difficult  enough,  but  there  are  other  complicating  factors. 

Mesoscale  atmospheric  dispersion  is  more  complicated  than  smaller-scale 
dispersion  because  the  mean  wind  field  cannot  be  considered  steady  or  horizontally 


9 


homogeneous  over  the  time  and  space  scales.  Wind  shear  plays  a  much  greater  role: 
horizontal  dispersion  can  be  enhanced  and  often  dominated  by  vertical  wind  shear  on 
smaller  scales  through  interaction  of  horizontal  differential  advection  and  vertical  mixing 
(Moran  and  Pielke,  1994:96).  Effects  of  mesoscale  phenomena  should  be  taken  into 
account  in  dispersion  studies,  but  this  cannot  be  done  without  complete  wind  and  rain 
field  data  (Rantalainen,  1993:143). 

Over  irregular  terrain  the  diffusion  and  transport  of  pollutants  are  significantly 
more  difficult  to  model  than  over  flat  terrain.  The  effect  of  the  topography  enhances 
diffusion  because  of  increased  mechanical  turbulence  and  wind  shear,  and  the  possibility 
of  a  larger  number  of  pollutant  trajectories  which  are  convoluted  by  topographic 
channeling  and  by  thermally  induced  flows  (Mullen  et  ah,  1978:188-191). 

Land-water  interfaces  mentioned  in  the  previous  section  may  be  further 
complicated  by  steep  planetary  boundary  layer  (PEL)  heights— some  on  the  order  of 
1500  m  over  a  grid  interval  of  10  km-associated  with  sea-breeze  fronts  and  enhanced  by 
the  topography.  In  mountainous  regions  with  large  scale  flows  a  significant  horizontal 
shift  in  the  maximum  PEL  height  relative  to  the  mountains  is  induced  by  a  corresponding 
displacement  of  the  thermal  ridge  (Lieman  and  Alpert,  1993:129). 

2.4  Example  of  Models 

In  the  course  of  generating  models  to  predict  weather  and  dispersion  in  the 
atmosphere,  many  different  techniques  have  been  generated.  It  is  interesting  to  compare 
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some  of  these  methods  and  see  the  limitations  and  benefits  derived  from  them  with 
OMEGA. 

The  traditional  source-oriented  approach  consists  of  solving  model  equations 
forward  in  time  for  given  emission  sources  of  pollutant  to  obtain  a  time  and  space 
distributed  concentration  field.  It  allows  calculation  of  various  air  pollution 
characteristics  for  any  number  of  receptors.  For  any  new  emission  scenario,  however,  the 
model  solution  must  be  repeated  (Uliasz  et  al.,  1994:104). 

With  any  grid,  available  data  does  not  lie  directly  on  the  grid  points.  A  common 
method  is  to  interpolate  to  the  model  grid  at  each  observation  time  by  an  optimal 
interpolation  method.  The  spatial  weighting  function  is  based  on  the  estimation  variance 
of  the  interpolated  data  obtained  from  this  method.  If  an  observation  is  located  precisely 
at  a  model  grid  point,  W  =  1,  and  if  a  model  point  is  sufficiently  distant  from  the 
observation  location,  W  =  0;  otherwise,  the  weight  varies  depending  on  spatial 
distribution  and  number  of  observations. 

The  following  subsections  under  2.4  describe  eight  different  models.  The  Newton 
relaxation  method  describes  a  means  of  predicting  flow  fields  using  limited  amounts  of 
data.  The  Gaussian  Plume  Model  is  described  since  it  is  generally  used  by  the  U.S.  EPA. 
The  Air  Resources  Laboratory  (ARL)  model  demonstrates  means  of  predicting 
deposition.  The  remaining  models  were  used  in  the  ATMES  exercise  and  presented  for 
comparison. 
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2.4. 1  Newton  Relaxation 


Newton  relaxation  used  to  predict  flow  fields  have  two  major  advantages.  One  is 
its  conceptual  simplicity  and  low  computational  demand.  Another  is  it  produces  a  flow 
field  that  is  a  combination  of  the  predicted  and  observed  variables  when  and  where  the 
observations  may  occur.  In  data  sparse  regions,  only  the  model  governing  equations  are 
used  to  predict  the  flow  field.  In  theory,  continuous  four-dimensional  data  assimilation 
should  be  superior  to  mass-consistent  diagnostic  models  that  simply  interpolate  the  wind 
field  in  three  dimensions  into  data  sparse  regions.  Spatial  and  temporal  interpolation 
techniques  employed  by  these  models  may  result  in  significant  wind  speed  and  direction 
errors,  especially  in  highly  complex  terrain. 

The  disadvantages  of  Newtonian  relaxation  are:  1)  the  coefficient  of  relaxation  is 
a  guess;  and  2)  assimilation  of  local  or  unrepresentative  components  can  occur  (i.e., 
microscale  observations  may  be  spread  over  a  relatively  large  area).  For  instance,  local 
influence  of  a  valley  floor  observation  can  be  spread  up  the  valley  walls  (Fast  and 
O’Steen,  1994:310-311). 

2.4.2  Gaussian  Plume  Models 

For  a  Gaussian  emission  model  assuming  continuous  emission,  mass 
conservation,  and  steady-state,  the  crosswind/vertical  concentration  distance  is  given  by 
using  the  standard  deviations  of  plume  distribution  for  all  directions:  Ox,  Gy,  and  Gz. 
These  values  depend  on  the  weather  parameters  and  stability  classes.  The  Gaussian 
distribution  is  given  by  the  ordinate  value  as  shown  in  Figure  1. 
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Figure  1 .  Gaussian  Distribution  Curve. 

(Turner,  1994:2-1,2-4). 

In  its  simplest  form,  the  Gaussian  plume  equation  for  downwind  concentrations  at 
ground  level  is  written; 

X(x,y,z=0)  =  Q  [KuOyOz]  *^^exp[(yVa^y  +  H  /a  z)]/2, 
where  x(x,y,z=0)  is  ground-level  concentration  at  downwind  distance  x  and  crosswind 
distance  y,  u  =  mean  wind  speed,  Q  =  source  intensity,  H  =  effective  height  of  release, 
and  Oy  and  Oz  are  lateral  and  vertical  coefficients  of  dispersion.  Plume  rise  and 
coefficients  of  dispersion  can  seldom  be  measured,  which  forces  one  to  resort  to  the  use 
of  empirical  formulae,  of  which  there  are  many  proposed.  The  variability  inherent  in 
these  parameters  warrant  careful  consideration  of  any  Gaussian  model,  and  it  is  perhaps 
advisable  to  carry  out  runs  using  several  of  the  best-performing  sigma  sets  (e.g.  the 
Brigg’s  Formulae  perhaps  most  widely  used)  to  get  a  feel  for  the  range  of  variability  to  be 
expected  (Carrascal,  1993:147,  156).  In  a  report,  Pasquill  reminded  the  U.S.  EPA  to 
make  the  best  use  of  Gaussian  modeling  through  specific  measurements  relating  closely 
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to  dispersion.  He  reemphasized  the  use  of  standard  deviation  of  the  wind  direction  angle 
determined  over  1-hr  periods. 

Accuracy  of  Gaussian  plume  models  are  limited  by  several  factors.  Errors  in  the 
emission  rate  will  propagate  directly  into  an  error  in  the  calculated  concentration.  Since 
wind  speed  generally  increases  with  height  above  the  ground,  estimation  of  the  wind 
speed  at  the  point  of  release  may  be  in  error  on  the  order  of  10  to  15  percent.  Since  the 
effective  height  of  release  is  dependent  upon  wind  speed  and  stack  parameters,  it  may 
also  have  errors  of  up  to  20  to  25  percent.  Dispersion  parameters  found  by  Pasquill 
stabilities  and  Pasquill-Gifford  dispersion  parameters  consider  the  atmosphere  in  six 
classes,  while  in  reality  it  is  a  continuum.  Considerably  different  concentrations  are 
calculated  with  a  change  of  one  stability  class  in  the  assumptions.  Larger  errors  are 
expected  with  extremes  of  stability  and  large  distances. 

At  the  source  of  release,  slight  errors  in  the  estimation  of  wind  direction  can  result 
in  tremendous  errors  of  concentration  where  the  problem  is  to  estimate  the  concentration 
at  specific  locations.  Often  in  this  case  the  magnitude  of  the  highest  downwind 
concentrations  under  stated  stability  and  wind  speed  are  estimated  quite  well,  but  the 
location  may  be  in  error.  Therefore,  it  is  important  to  make  exceptionally  good  estimates 
of  the  wind  direction  for  each  time  period  or  expect  to  put  up  with  large  error  bounds, 
perhaps  as  much  as  a  factor  of  ten. 

Gaussian  models  are  good  for  estimation,  but  their  application  is  limited  by  the 
assumptions  not  applying  to  many  release  scenarios  and  by  the  accuracy  in  measuring  the 
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model  parameters.  Its  simplicity  and  ability  to  predict  many  short  term  velocities  by  a 
normal  distribution  are  still  attractive  (Turner,  1994:2-13,  2-15,  5-1). 


2.4.3  Air  Resources  Laboratories  (ARL) 

ARL  developed  a  computerized  model  to  calculate  transport,  diffusion,  and 
deposition  of  effluents  on  regional  and  continental  scales.  Each  puff  diffuses  according 
to: 

Cm  =  (Q/2pah^Zm)exp(-R^/2ah^) 

where  Cm  is  ambient  air  concentration  in  the  mixed  layer,  Q  is  emission  amount  per  puff. 
Oh  is  horizontal  standard  deviation,  Zm  is  height  of  the  mixed  layer,  and  R  is  distance 
from  puff  center.  The  concept  of  a  deposition  velocity  is  used  to  calculate  dry  deposition 
along  a  trajectory  and  an  empirical  scavenging  ratio  is  used  for  wet  deposition.  The 
fraction  of  mass  removed  from  the  mixed  layer  by  dry  deposition  is: 

CmVddt/CmZm 

where  Vj  is  the  dry  deposition  velocity  and  dt  is  the  time  interval  at  which  puff 
concentrations  are  calculated  (Heffter  and  Ferber,  1990:400-407). 

2.4.4  WSPEEDI 

The  Worldwide  version  of  System  for  Prediction  of  Environmental  Emergency 
Dose  Information  (WSPEEDI)  was  revised  to  include  terrain  following  vertical 
coordinates  with  compressibility  of  the  atmosphere  being  considered.  A  particle 
dispersion  model  is  used  to  simulate  the  atmospheric  transport  of  radionuclides.  The 
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radioactive  plume  is  expressed  by  a  mass  of  particles  with  separate  calculations  for  the 
position  of  each  particle.  The  equations  incorporate  effects  of  eddy  diffusion  and 
horizontal  diffusion  to  describe  particle  movements. 

WSPEEDI  was  a  participant  in  the  ATMES  exercise.  The  computed  surface  air 
concentration  of  Cs-137  was  directly  compared  with  measurements  at  18  sites. 

Agreement  is  generally  good  in  Central  Europe  and  Kozanice  (Greece).  At  some  sites  in 
Western  Europe,  the  simulated  results  correspond  to  the  measurement  well  during  several 
days  after  the  arrival,  but  the  computed  concentration  decreases  rapidly  after  that.  This 
rapid  decrease  was  due  to  the  westerly  winds  which  prevailed  over  these  areas  bringing  in 
clean  air  after  May  3.  This  suggests  the  possibility  the  plume  reached  further  west  than 
simulated.  The  uncertainty  reveals  the  difficulty  in  pinpointing  values  of  deposition, 
which  greatly  depend  on  local  precipitation  (Ishikawa,  1994:969-978). 

2.4.5  Canadian  Tracer  Model 

Extensive  development  work  with  a  three-dimensional  (3-D)  tracer  model  was 
conducted  by  the  Canadian  Meteorological  Center.  The  new  model  was  tested  on  the 
Chernobyl  case  using  objectively  analyzed  meteorological  fields  over  a  period  of  one 
month.  The  major  limitation  of  the  system  used  in  these  studies  was  the  use  of  observed 
meteorological  fields.  Thus,  the  development  of  a  predictive  version  of  the  tracer  model 
is  a  natural  continuation  of  previous  work. 

The  models  which  simulate  atmospheric  tracers  are  relatively  complex  because  of 
the  need  to  represent  a  wide  range  of  physical  and  chemical  processes.  To  simulate  the 
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interaction  between  meteorological  processes  and  atmospheric  transport  of 
nonconservative  species,  one  needs  a  set  of  thermodynamic  equations  coupled  with  the 
mass  conservation  equations  for  atmospheric  tracers.  The  system  in  the  most  general 
form  is  as  follows: 


dU/dt  =  F  -  (1/p)  Vp  -  nu 

(1) 

d  In  p/dt  =  -  V  •  U 

(2) 

dT/dt  =  At  (L/cp)  Q  -  (L/cp)Eo 

(3) 

p  =  pRT 

(4) 

dq/dt  =  -Q  -1-  Eo 

(5) 

dm/dt  =  Q  -  (P  -  Er)  -  Eo 

(6) 

dnVdt  =  S‘  (n', ...,  m,  T) 

(1=1,...,  M) 

(7) 

Where  the  symbols  have  the  following  meanings: 

U- velocity  field,  F-mass  forces,  p-atmospheric  pressure,  0  -  operator  representing 
turbulent  mixing,  p  -density,  V  -  gradient  operator,  T-temperature,  R-gas  constant, 
q  -  mixing  ratio  of  water  vapor,  L  -  latent  heat  constant.  At  -  dynaimc  part  of  tendency 
rate  for  temperature,  Q  -  latent  heat  of  condensation,  Eo  -  evaporation,  m  -  cloud  water, 
P  -  precipitation  rate,  Er  -  evaporation  of  rain,  n‘  -  mixing  ratio  for  atmospheric  trace 
species,  S'( )  -  source  and  sinks  for  atmospheric  trace  species,  d/dt  -  material  derivative 
(d/dt  =  3/3t  +  U  •  V),  M  -  number  of  atmospheric  trace  species  simulated  in  the  model. 
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The  system  (l)-(7)  consists  of  (8  +  M)  equations  and  simulates  the  meteorological 
processes  and  chemical  reactions  of  atmospheric  trace  species.  The  boundary  conditions 
for  the  dynamic  and  chemical  part  of  the  system  have  to  be  specified,  and  they  depend  on 
the  particular  geometry  of  the  domain  being  considered. 

The  first  eight  equations  of  the  system  are  in  fact  a  non-hydrostatic  meteorological 
model.  The  model  equations  include  the  bulk  representation  of  microphysical  processes 
such  as  condensation,  the  formation  of  clouds  and  precipitation  and  the  evaporation  of 
rain.  The  atmospheric  boundary  layer  processes  are  represented  by  the  operator  0  which 
takes  into  account  the  momentum  transfer  by  convective  processes  and  by  dynamically 
generated  turbulence.  The  set  of  Eq.  (7)  represents  a  nonlinear  atmospheric  chemistry. 
The  set  given  by  Eqs.  (l)-(7)  is  designed  mostly  for  application  on  a  regional  scale  and  is 
far  too  complex  for  most  applications  involving  the  simulation  of  atmospheric  processes 
in  a  hemispheric  or  global  scale.  One  common  simplification  is  to  assume  hydrostatic 
equilibrium  in  a  tracer  model  driven  by  a  primitive  equations  meteorological  model. 

These  simplifications  affect  calculations  of  turbulent  vertical  fluxes  of 
momentum,  heat  and  moisture  in  the  surface  layer,  and  scale  horizontal  mixing.  Stable 
precipitation  is  currently  calculated  by  removing  the  excess  of  moisture  above  some 
threshold  value  of  relative  humidity.  The  moist  convection  is  parameterized.  Cloud 
cover  required  in  the  calculation  of  radiative  heating  is  obtained  as  a  function  of  the  local 
relative  humidity.  Surface  temperature  in  the  spectral  model  is  calculated  by  solving  an 
energy  balance  equation  over  land;  over  water,  surface  temperature  is  obtained  from 
climatological  data.  Wet  scavenging  is  represented  by  a  statistical  parameterization. 
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The  source  term  for  Chernobyl,  Q‘ ,  was  approximated  using  the  concept  of  the  3- 
dimensional  virtual  source.  The  release  in  this  approximation  is  represented  by  a  function 
of  three  spatial  variables  and  time: 

Q‘  (tlx,  riy,  a,  t)  =  ((E'(t)  F(a))/27tc'H))Exp[-rWH]  (8) 

where: 

E'(t)  —  amount  of  radioactivity  released 

-1- 1 

■<Jb 

go 

<5  Y 

f(a)  -  function  describing  the  vertical  distribution  of  the  release 
(t|°x,  T|°y)  -  coordinates  of  the  source  of  the  release 
r  =  ((T|°x -(rix)^  +  (Tl°y -Tly)^)’'^ 

Oh  -  standard  deviation  of  the  horizontal  mass  distribution 
as,  ctt,  are  values  at  the  bottom  and  top  of  the  domain,  respectively 
R  -  specific  gas  constant  of  air 
T  -  temperature 

The  effective  release  was  a  function  of  time  with  height  changing  from  4000  m 
just  after  the  explosion  to  1500  m  subsequently.  The  initial  height  of  4000  m  was 
assumed  to  be  in  direct  response  to  the  initial  explosion.  The  smaller  vertical  extension 
of  1500  m  which  was  assumed  during  the  days  following  the  accident  mostly  reflects  the 
intensive  convection  over  the  burning  reactor  and  subgrid-scale  vertical  mixing  of 
radioactive  material  in  the  lower  part  of  the  troposphere. 


F(a)  -f(a)-|— |- 
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The  tracer  was  solved  initially  using  data  stored  in  the  archiving  system  of  the 
Objective  Analysis  at  the  Canadian  Meteorological  Center.  The  analysis  of 
meteorological  fields  at  the  time  of  the  accident  was  performed  every  6  hours  on  a 
Gaussian  latitude-longitude  grid  with  a  resolution  of  128*32  points  over  the  Northern 
Hemisphere.  The  meteorological  fields  for  intermediate  time  levels  were  obtained  by 
time  interpolation. 

The  verification  of  the  results  from  the  diagnostic  model  was  performed  for  a 
network  covering  the  Northern  Hemisphere  and  indicates  that  the  simulation  was  quite 
accurate.  Predicted  times  of  arrival  and  the  times  of  arrival  of  maximum  activity  agree 
quite  well  with  the  observed  values.  The  ratio  of  the  calculated  and  observed  values 
varies  between  0.21  and  2.80  with  a  mean  value  of  1.05. 

The  results  of  the  diagnostic  calculations  indicate  that  during  the  first  two  days 
following  the  Chernobyl  accident,  the  radioactive  cloud  was  transported  mostly  towards 
Scandinavia.  A  second  southern  segment  of  the  cloud  had  spread  southeastward  over  the 
Black  Sea  in  the  direction  of  the  Middle  East.  On  April  29,  a  well  developed  westerly 
flow  began  to  transport  radioactive  material  across  the  former  USSR. 

The  2  day  prediction  has  relatively  good  accuracy,  as  the  model  was  quite  good  in 
predicting  the  main  directions  of  the  transport.  The  quality  of  the  4  day  prediction, 
however,  deteriorated  substantially.  The  deficiency  of  the  predictive  model  is  manifested 
mostly  by  the  lack  of  radioactive  material  over  western  Europe.  The  first  run  of  the 
predictive  model  showed  good  prediction  on  a  regional  scale  for  the  first  2-3  days  after 
the  release  (Pudykiewicz,  1990:213-225). 
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Another  approach  to  Chernobyl  involved  replacing  the  usual  Pasquill  dispersion 
coefficients  (sigmas)  for  the  simplified  Gaussian  plume  diffusion  approach  in  this  terrain 
with  dispersion  coefficients  related  empirically  to  the  terrain  roughness  and  atmospheric 
turbulence.  A  steady  state  modeling  approach  is  used;  neglect  of  transport  time  simplifies 
the  model  considerably  but  causes  the  model  to  overreact  to  changes  in  meteorology 
when  modeling  receptors  are  long  distances  from  the  sources.  In  spite  of  this 
shortcoming,  the  model  predicts  reasonably  well  for  all  averaging  times  (Mullen  et  al., 
1978:188-191). 

2.4.6  MATHEW 

A  Lagrangian,  mass-adjusted,  three-dimensional  wind  field  model  (MATHEW) 
was  developed  to  provide  wind  fields  for  a  pollutant  transport  model  (ADPIC).  The  two 
basic  boundary  conditions  imposed  on  the  field  are  constant  mass  flux  and  zero  mass  flux 
corresponding  to  reflection  on  the  boundaries  (Klug,  1992:  68).  The  wind  model 
incorporates  terrain  explicitly,  is  site  independent,  uses  available  meteorological 
measurements,  is  stable,  and  calculates  three-component  velocities  at  a  large  number  of 
grid  points.  Constant  air  density  was  assumed,  but  variable  air  density  can  be  included 
with  little  modification. 

Comparison  of  results  with  field  data  taken  during  tracer  releases  shows  the  model 
calculations  to  be  within  a  factor  of  2  in  50%  of  the  comparisons  and  within  an  order  of 
magnitude  for  90%  of  the  tests.  Examination  of  model  calculations  shows  the  root-mean- 
square  errors  to  be  10%  in  wind  speed  and  5-10°  in  wind  direction,  except  near 


21 


topographic  barriers  where  the  errors  increase  to  follow  flow  around  and  over  the  terrain 
features  (Sherman,  1978:312-319). 

2.4.7  MEDIA 

This  model  involves  numerical  schemes  that  can  deal  with  advection  and 
diffusion  processes  without  distortion  by  numerical  effects.  Pollution  concentration  is 
advected  by  the  wind  field  as  predicted  by  the  operational  coupling  model.  To  simplify 
the  procedure,  the  diffusion  is  modeled  using  exchange  coefficients.  It  means  pollutant  is 
assumed  to  be  diffused  in  the  same  way  as  water  vapor. 

Sinks  for  the  model  are  treated  as  follows:  Wet  deposition  due  to  scavenging  by 
precipitation  is  computed  using  a  global  coefficient  of  air-to-water  transfer,  which 
roughly  describes  dilution  or  catching.  Dry  deposition  is  modeled  using  a  coefficient 
dimensionally  equal  to  a  deposition  velocity  multiplied  by  the  concentration  in  air  near 
the  ground.  Radioactive  decay  is  another  sink.  Diffusion  is  modeled  on  a  subgrid  scale 
by  a  Gaussian  distribution.  On  the  six  sides  of  the  integration  area  (four  lateral  sides,  top 
of  the  atmosphere  and  ground  level)  boundary  conditions  are  treated  as  outgoing  fluxes 
according  to  Orlanski  (1976). 

The  data  used  was  from  the  European  Centre  for  Medium  Range  Weather 
Forecasts.  Time-evolution  was  according  to  a  step  function  as  described  in  the  Technical 
Specifications  Document  of  the  ATMES  Report. 
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Results  were  very  satisfactory  for  all  but  one  location.  Problems  with  monitoring 
stations  close  to  the  accident  may  be  the  result  of  linear  interpolation  over  6  hours  of 
large-scale  meteorological  fields  or  errors  in  knowledge  of  the  source  term. 

Sensitivity  tests  were  performed  on  the  model.  A  space  resolution  of  2  x  2 
degrees  did  not  result  in  a  great  loss  of  quality  on  the  average  values.  Reducing  the  time 
step  to  1  hour  instead  of  3  hours  did  not  notably  improve  the  results.  Reducing  horizontal 
diffusion  leads  to  a  great  loss  of  quality  on  the  smaller  values  without  any  benefit  for  the 
extrema^  increasing  it  leads  to  smooth  minima  and  maxima.  The  chosen  value  seemed 
best. 

Sensitivity  to  deposition  velocity  is  rather  important  and  shows  it  is  necessary  to 
know  and  accurately  describe  the  characteristics  of  the  transported  material  and 
underlying  surface  over  which  it  will  be  deposited.  Omitting  rainfall  gives  almost  the 
same  results  for  ambient  concentration.  This  can  be  explained  by  the  fact  that  during  the 
10  days  after  the  accident,  the  radioactive  cloud  was  seldom  affected  by  washout. 

Average  concentrations  in  air  are  not  the  best  indicators  since  the  average  losses  due  to 
washout  is  low  compared  to  the  amount  of  pollutant  released  during  such  a  short  time 
(Piedelievre  et  al.,  1990:1205-1220). 


2.4.8  LORAN 

LORAN  (Long  range  atmospheric  advection  of  nuclides)  makes  use  of  horizontal 
wind  field  varying  linearly  in  time  and  space.  The  advection  velocity  for  the  current  time 
interval  dt  at  point  (x,y)  for  the  given  geopotential  height  used  is  calculated  via  linear 
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interpolation  of  the  comer  point  values  and  for  the  times  tl  and  t2.  A  typical  value  of  dt 
is  1  hr. 

Dry  and  wet  deposition  is  evaluated  by  means  of  an  average  deposition  value  for 
each  pollutant  or  nuclide  used  in  the  simulation.  In  principle,  deposition  velocity  should 
depend  on  atmospheric  turbulence  and  surface  resistance.  The  long  duration  of  the 
release,  however,  and  the  lack  of  information  on  surface  resistance,  has  prompted  the  use 
of  the  value  of  0.003  m/s.  Scavenging  of  contaminated  air  by  rain  is  described  by  a 
scavenging  coefficient  times  the  concentration. 

Results  of  the  comparison  for  air  show  a  constant  overestimation  of  the 
concentration.  This  can  be  explained  by  the  fact  that  a  considerable  fraction  of  the  source 
did  not  travel  to  long  distances  from  the  release  point.  Results  for  cumulative  deposition 
show  good  agreement  in  the  two  cumulative  distributions  of  results.  The  overall  bias  for 
deposition  again  indicates  overestimation. 

The  model  LORAN  would  have  been  among  the  best  had  it  participated  in  the 
ATMES  exercise.  It  is  believed  that  its  good  performance  is  mainly  due  to  its 
prescription  of  the  mixing  layer  growth.  Results  of  the  comparison  of  cumulative 
deposition  data  are  very  different  when  data  close  to  the  source  point  are  included.  This 
indicates  that  the  model  does  not  behave  well  close  to  the  source;  probably  due  to  the 
simple  assumption  on  the  constant  reference  level  and  on  instantaneous  air  mixing  in  this 
level.  The  model  actually  seems  to  perform  better  if  Russian  deposition  data  is  not 
included  (Galmarini  et  al.,  1992:143-154). 
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2.5  Evaluation  of  Models 


Model  accuracy  is  difficult  to  determine  since  it  depends  upon  both  the  ability  of 
the  model  to  simulate  the  physics  of  the  atmosphere  and  the  accuracy  of  the  input 
information.  Most  evaluations  contain  the  product  of  these  two  effects.  Although  a 
comparison  between  model  estimates  and  measurements  may  verify  a  model,  a  good 
result  may  be  due  to  compensating  errors  in  various  portions  of  the  model.  Thus, 
independent  verification  of  portions  of  the  model  is  highly  desirable.  The  comparison 
may  be  subjective  by  displays  of  scatter  plots  or  quantitative  such  as  deriving  various 
statistics  from  data  sets.  The  use  of  laboratory  data  and  numerical  simulations  of 
dispersion  is  stressed  because  they  can  be  repeated  to  generate  well-defined  ensemble- 
mean  concentration  fields. 

Two  questions  of  plume  dispersion  models  are  of  primary  concern:  1)  How  well 
does  a  model  predict  the  high  ground-level  concentrations;  2)  Is  the  model  based  on 
sound  physical  principles  and  give  good  predictions  for  the  “right”  reasons? 

There  are  two  major  steps  in  a  model  physical  evaluation.  The  first  is  an 
assessment  of  the  scientific  formulation  and  modeling  assumptions  for  the  physical 
problem  of  interest.  The  second  is  an  evaluation  of  model  predictions  using 
measurements  from  laboratory  and  field  experiments. 

Many  different  forms  of  tests  are  used  for  model  evaluation.  Performance 
evaluation  tests  the  performance  under  different  conditions  through  a  partitioning  of  the 
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data.  Sensitivity  testing  determines  the  changes  in  model  output  due  to  specific  changes 
in  an  input  parameter  (Turner,  1979:502-519). 

A  key  step  in  evaluation  is  separating  the  model  error— the  difference  between  the 
observed  and  predicted  values  of  concentration-from  the  variability  and  other  sources  of 
concentration  variance.  This  is  best  done  through  a  residual  analysis  in  which  one 
examines  the  statistics  of  the  residual  between  the  observed  and  predicted  concentrations. 
The  key  objective  in  residual  analysis  is  to  show  the  difference  exhibits  no  trends  with 
any  of  the  ensemble  variables  (Weil,  1994:224). 

The  following  statistical  tests  are  used  for  comparison  of  model  results  with  actual 
data:  True  difference,  normalized  true  difference,  absolute  difference,  absolute  fractional 
bias,  variance  of  true  difference,  variance  of  fractional  bias,  mean  square  error,  and 
correlation  coefficient  (Ciolek,  1994:237-238). 

Model  evaluation  for  the  Chernobyl  Nuclear  Accident  is  based  on  the  comparison 
between  the  following  observed  and  computed  analysis:  Cs-137  air  concentration 
samplings  paired  in  space  and  time;  time-integrated  concentrations  at  each  locality;  time 
of  arrival  of  the  cloud.  The  model  simulation  and  consequently  all  the  analysis  are 
limited  to  the  first  12  days  after  the  beginning  of  the  release,  i.e.,  until  24:00  GMT  on 
7  May. 

The  APOLLOS  model  defined  the  computed  arrival  time  of  the  cloud  at  a  certain 
location  as  the  beginning  of  the  first  time  step  at  which  the  Cs-137  concentration  exceeds 
a  threshold  value  CT;  the  observed  arrival  time  is  defined  as  the  initial  time  of  the 
sampling  interval  with  measured  concentration  exceeding  CT.  A  threshold  value 


26 


CT  =  0.1  Bq  m'^  was  chosen,  with  the  exception  of  a  few  sites  whose  first  measurement 
was  higher  than  0.1,  for  which  CT  =  1  Bq  Due  to  large  uncertainty  of  the  observed 
arrival  time  and  the  unavoidable  arbitrary  and  inhomogeneous  choice  of  the  threshold, 
this  part  of  the  study  is  intended  to  give  only  a  rough  indication  of  the  model  capabilities 
of  estimating  the  time  of  transport  of  the  cloud. 

There  are  several  potential  reasons  for  inaccuracy  of  model  results  when  a  long- 
range  and  long-duration  dispersion  episode  like  the  Chernobyl  release  is  simulated. 

Among  others,  the  most  important  are  perhaps  the  prescription  for  source  rate  and  shape, 
the  three-dimensional  wind  field  used  for  calculating  the  particle  trajectories,  the 
parameterizations  of  boundary-layer  height  and  the  horizontal  diffusion,  the  treatment  of 
the  cloud  scavenging  by  dry  and  wet  deposition.  Besides,  the  relative  weight  of  each  of 
these  sources  of  uncertainty  strongly  depends  on  the  method  used  for  comparing  the 
results. 

After  several  air  quality  model  evaluation  exercises  involving  a  large  number  of 
source  scenarios  and  types  of  models,  it  is  clear  that  the  magnitudes  of  the  uncertainties  in 
model  predictions  are  similar  from  one  application  to  another.  For  example,  when 
considering  continuous  point  sources  and  receptors  at  distances  of  about  0. 1  km  to  1  km 
downwind,  uncertainties  in  ground-level  concentration  predictions  lead  to  typical  mean 
biases  of  about  ±20  to  40%  and  typical  relative  root-mean-square  errors  of  about  60  to 
80%.  It  is  not  uncommon  to  see  overprediction  by  50%  at  one  site  and  underprediction  of 
50%  at  a  second  site  for  two  otherwise  identical  model  applications.  This  fundamental 
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level  of  model  uncertainty  is  likely  to  exist  due  to  data  input  errors  and  stochastic 
fluctuations,  no  matter  how  sophisticated  a  model  becomes  (Hanna,  1993:3). 
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3.  Methodology 


3. 1  Use  of  Medium  Range  Forecast  Data  over  the  U.S. 

Running  OMEGA  using  a  domain  within  the  U.S.  will  provide  means  of  verifying 
its  ability  to  predict  wind  fields.  We  selected  the  Medium  Range  Forecast  (MRF)  data 
from  the  National  Weather  Service  (OMEGA  data  processor  also  comes  with  the 
capability  to  use  Nested  Grid  Model  output  and  Naval  Oceanographic  Global 
Atmospheric  Prediction  System  (NOGAPS)  data  from  Fleet  Numerical  and 
Oceanographic  Center,  or  hand  entered  data).  The  data  includes  actual  measurements  at 
OOOOZ  of  the  day  of  interest  with  predictions  for  1200Z  and  2400Z.  The  ability  of 
OMEGA  to  predict  wind  fields  for  the  domain  of  interest  will  depend,  in  part,  on  the 
accuracy  of  the  later  time  MRF  predictions. 

The  domain  of  interest  spans  36.5°  N  by  43.5°  N  and  80.25°  W  to  89.75°  W  (see 
Figure  2  on  the  next  page  for  a  view  of  the  domain).  These  ranges  were  chosen  to  make 
best  use  of  the  10  km  resolution  terrain  data  included  with  the  model  while  presenting 
reasonable  run  times  for  the  workstation.  The  days  covered  the  16th  through  the  18th  of 
September,  1995.  These  were  chosen  because  of  a  trough  seen  in  the  wind  flow.  The 
boundary  conditions  will  then  have  wind  field  in  basically  the  same  direction,  but  the 
domain  (with  wind  fields  generated  by  OMEGA)  will  not  yet  account  for  the  trough.  The 
ability  of  OMEGA  to  generate  the  trough  within  the  boundary  will  partially  validate  its 
ability  to  predict  wind  fields. 
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Figure  2.  Domain  of  MRF  Data  Run. 
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OMEGA  is  run  three  separate  times  with  one  day  of  MRP  data  per  run.  This 
involves  analysis  data  at  OOOOZ  and  forecast  data  at  1200Z  and  2359Z.  Later,  the  run  may 
be  repeated  using  the  boundary  conditions  for  each  of  the  cases  in  sequence  over  a  longer 
period.  This  makes  use  of  analysis  data  alternating  with  forecast  data  every  twelve  hours. 
We  expect  the  latter  run  to  be  more  consistent  with  analysis  data  because  more  of  the 
boundary  conditions  use  analysis  data.  The  general  flow  of  the  actual  weather  patterns 
and  the  velocity  will  be  compared  with  OMEGA  predictions  at  four  levels.  Assessment 
at  the  300,  500,  700,  and  850  mb  levels  at  each  station  at  twelve-hour  intervals  leads  to 
80  data  readings.  OMEGA  predictions  involve  interpolation  because  the  layers  do  not 
exactly  correspond  to  these  pressures.  Therefore,  a  third-order  Lagrangian  polynomial 
using  two  data  points  on  either  side  of  the  value  desired  was  selected  to  obtain  the  values 
for  wind  speeds. 

The  means  of  assessing  the  predictions  of  OMEGA  comes  through  a  two-tailed 
t-test.  The  weather  stations  in  the  domain  available  for  us  from  the  National  Weather 
Service  analysis  maps  are  Dayton,  OH;  Huntington,  WV;  and  Pittsburgh,  PA.  The 
resolution  of  the  winds  from  the  analysis  maps  is  2.5  knots  and  the  resolution  of  the 
direction  is  about  10°.  The  u  and  v  components  estimated  from  the  analysis  map  vectors 
are  then  compared  with  the  predictions  of  OMEGA.  The  test  chosen  reflects  use  of  paired 
data  and  the  null  hypothesis  chosen  was  that  the  mean  of  the  wind  speeds  was  equal  for 
the  predicted  wind  and  actual  winds.  A  95%  confidence  interval  was  chosen  with  the  null 
hypothesis  rejected  when  the  limits  set  by  the  t-value  was  exceeded. 
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3.2  Verification  with  HIRAS  Data  and  the  Chernobyl  Nuclear  Accident 

3.2. 1  Background  from  HIRAS  USAFETAC  Climatic  Database  Users  Handbook 
No.  5  (October  1988:3,  8) 

Assessment  of  OMEGA  will  also  come  through  verification  with  Chernobyl  data. 
The  model  of  choice  is  HIRAS  data,  or  High  Resolution  Analysis  System,  generated  and 
obtained  from  archives  at  OL-A,  Air  Force  Combat  Climatology  Center  (AFCCC).  The 
data  is  recorded  at  six-hour  intervals  and  involves  only  analysis  data.  A  disadvantage, 
however,  is  the  limitations  in  the  1986  data  as  described  later  in  this  section.  Use  of  all 
analysis  data  would  lend  more  credibility  to  the  ability  of  OMEGA  to  predict  wind  fields 
based  on  observed  data  rather  than  forecast  data. 

OMEGA  does  not  support  use  of  archived  HIRAS  data  (the  newer  data  is  in  a 
different  format  and  would  be  used  in  future  scenarios),  thus  requiring  a  special  data 
preprocessor.  HIRAS  produces  global  upper-air  analyses  on  a  2.5°  by  2.5° 
latitude/longitude  grid.  All  standard  pressure  levels  are  available.  HIRAS  uses  a  variety 
of  observations  taken  from  land  stations,  ships,  buoys,  aircraft,  RAOBs,  PIBALs,  and 
satellites.  HIRAS  analyzes  five  elements— heights,  u-  and  v-component  winds, 
temperature,  and  relative  humidity-directly;  all  other  HIRAS  elements  are  derived  from 
these  five. 

HIRAS  has  two  main  components:  a  “first-guess”  and  an  analysis  model.  HIRAS 
uses  an  analysis  technique  called  “optimum  interpolation”,  or  01.  This  model  is  an 
adaptation  of  the  analysis  used  at  the  National  Meteorological  Center.  01  takes  into 
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account  three  factors;  1)  distance  between  observations  and  grid  points;  2)  accuracy  of 
the  observing  instruments;  and  3)  expected  accuracy  of  the  first-guess  value. 

Distance  between  observations  and  grid  points  is  incorporated  by  assigning 
weights  to  observations  surrounding  each  grid  point.  In  HIRAS,  these  weights  decrease 
exponentially  with  distance.  Each  observation  is  allowed  to  affect  the  first-guess  analysis 
depending  on  how  close  it  is  to  the  grid  point.  If  a  grid  point  has  observations  nearby,  the 
first-guess  value  is  corrected;  otherwise,  it  remains  unchanged. 

The  second  factor  manifests  a  major  advantage  of  01.  It  assigns  every  instrument 
type  a  unique  “expected  error”  that  has  been  determined  statistically.  The  basic  01  rule  is 
that  the  lower  the  expected  instrument  error,  the  more  weight  that  the  observation  will 
receive. 

The  last  factor  is  accuracy  of  the  first  guess.  Each  HIRAS  analysis  produces  two 
types  of  fields:  analyses  and  errors.  Analyses  are  the  standard  grid-point  analyses 
previously  discussed.  Error  fields,  however,  are  specialized  fields  unique  to  01.  The 
error  fields  are  used  as  “running”  standard  deviations,  and  indicate  how  accurate  the 
analysis  is  at  each  grid  point.  The  more  observations  available,  the  better  the  analysis  and 
the  lower  the  expected  error.  Over  data-rich  areas  like  North  America  and  Europe,  error 
values  are  low. 

HIRAS  upper-air  analysis  provides  all  the  initial  conditions  for  AFCCC  models, 
especially  the  Global  Spectral  Model.  The  main  purpose  of  an  initialization  scheme  is  to 
control  imbalances  in  initial  conditions  which  cause  development  of  motions  (like  gravity 
waves)  that  the  model  cannot  resolve.  These  perturbations  could  be  real  (from  micro-  or 
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sub-grid  scale  phenomena)  or  fictitious  (from  small  observational  errors).  Nearly  every 
numerical  forecast  model  uses  some  form  of  initialization  to  compensate  for  these. 

Quality  control  applied  at  AFCCC  includes  manual  and  automated  checks. 
Forecasters  perform  the  manual  checks  by  adding  bogus  observations  or  deleting  real 
observations.  If  they  believe  anything  to  be  erroneous  while  reviewing  the  first-guess 
model  output,  forecasters  can  introduce  artificial  data  for  use  in  subsequent  analysis.  For 
example,  if  they  think  the  first-guess  has  moved  an  upper-air  trough  too  far  east,  artificial 
data  may  be  introduced  to  move  it  back  to  where  it  belongs. 

Automated  quality  control  is  directed  at  throwing  out  bad  observations.  The  first 
steps  are  when  HIRAS  compares  each  observation  with  the  first-guess  analysis.  If  an 
observation  grossly  disagrees  with  the  first-guess  (defined  to  be  when  the  difference 
between  the  observed  value  minus  the  first  guess  value  exceeds  the  analysis  error  field  for 
the  analysis  point  under  consideration  by  more  than  eight  error  standard  deviations),  it  is 
immediately  rejected.  If  an  observation  just  barely  passes  the  gross  checks,  it  is  flagged 
and  submitted  to  a  second  procedure  that  compares  it  with  a  nearby  observation.  If  one  of 
these  flagged  observations  significantly  disagrees  with  its  neighboring  observation 
(defined  to  be  when  the  difference  between  the  two  values  is  more  than  four  times  the 
analysis  error  at  the  analysis  point  being  considered),  it  is  rejected. 

There  are  some  known  problems  with  data  that  exist  for  the  Chernobyl  case. 
Moisture  analyses  are  not  from  HIRAS,  but  from  the  old  MULTAN  model.  Data 
elements  vorticity,  precipitable  water,  dew  point  depression,  and  relative  humidity  are  not 
available.  This  should  not  affect  the  performance  of  the  OMEGA  model  but  prevents  a 
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method  of  checking  the  model  performance.  Surface  temperature  is  actually  an 
extrapolation  from  the  1,000-miliibar  temperature  analysis.  This  data  is  not  surface 
temperature  and  should  not  be  used  as  such.  In  addition,  surface  elements  u-wind,  v- 
wind,  and  specific  humidity  are  not  available. 

3.2.2  Atmospheric  Transport  Model  Evaluation  Study  (ATMES)  Requirements 

Validation  using  Chernobyl  Nuclear  Accident  data  ideally  would  follow  the 

requirements  stated  in  the  Technical  Specification  Document  used  for  the  ATMES 
exercise.  The  exercise  encompassed  the  European  region  between  100°  W  to  40°  E  and 
35°  N  to  70°  N.  A  smaller  area  was  addressed  in  the  ATMES  exercise  by  limits  4°  E  to 
36°  E  and  43°  N  to  62°  N.  Model  results  were  requested  for  a  14  day  period,  with  model 
starting  time  reflecting  the  Chernobyl  source  term  data  viewed  as  a  step  function. 
Estimated  effective  height  of  the  initial  plume  center-of-mass  is  also  given.  Models 
predicted  several  parameters  including  surface  air  concentration,  arrival  times  for  Berlin 
and  Munich,  and  Cs-137  deposition  rates. 

3.2.3  Methods  Used  for  OMEGA 

The  domain  chosen  for  OMEGA  validation  will  allow  a  proof-of-concept.  The 
grid  chosen  was  10°  E  to  34°  E  and  47°  N  to  62°  N  (see  Figure  3  on  next  page).  This  was 
chosen  to  include  in  the  domain  most  of  the  area  of  deposition  found  in  the  first  four 
days.  The  resolution  chosen  was  30  to  80  km  using  the  stretchable  horizontal  grid.  This 
is  comparable  to  resolutions  used  by  ATMES  exercise  participants.  The  HIRAS  data 
used  for  the  model  was  from  25  April  1986,  OOZ,  to  29  April  1986,  OOZ.  Modeling  the 
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Figure  3.  Domain  of  HIRAS  Data  Run. 


source  term  required  writing  a  Fortran  90  program  according  to  the  ATMES  Technical 
Specifications  Document. 

Several  factors  required  a  limited  scope  for  the  Chernobyl  validation.  The  number 
of  cells  in  OMEGA  is  limited  to  around  6,000.  In  order  to  keep  the  resolution 
comparable  to  other  models,  the  domain  must  be  constrained.  The  existing  IBM  AIX 
system  runs  approximately  the  same  length  of  time  as  the  simulation  time  for  this  part  of 
the  validation.  Output  files  generated  every  hour  require  12  MB  of  space.  Also,  the  data 
files  run  through  the  preprocessor  took  great  disk  space.  These  constraints  limited  the 
area  of  domain  and  the  duration  of  the  simulation. 

Limitations  in  the  OMEGA  model  itself  would  not  easily  lend  it  to  the 
comparison  of  ATMES.  As  mentioned  in  the  first  chapter,  OMEGA  has  different  ways 
of  tracking  releases.  The  puffs  released  would  follow  Lagrangian  treatment  in  tracer 
mode  of  the  model.  This  would  track  the  centroid  of  the  puffs  but  require  extensive 
integration  to  determine  the  concentration  along  each  time  step  and  the  amount  deposited 
to  the  ground.  This  prevented  use  of  the  user  defined  source  term  as  intended,  but  further 
efforts  to  validate  OMEGA  may  make  use  of  the  code  (see  Appendix  D).  The  Eulerian 
tracer  capability  permits  release  of  a  specified  concentration  and  allows  it  to  disperse,  but 
this  mode  does  not  support  a  varying  release  source  required  to  simulate  Chernobyl. 
Therefore,  the  validation  is  limited  to  proof-of-concept  that  the  wind  fields  generated  and 
the  deposition  footprint  predicted  by  the  tracers  show  the  capability  to  adequately  portray 
a  release  such  as  Chernobyl. 
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The  analysis  of  the  wind  fields  is  similar  to  that  of  the  runs  over  the  U.S.  Values 
used,  though,  are  at  ground  level,  850  mb,  and  500  mb,  and  only  at  OOZ  at  each  day.  The 
following  cities  were  chosen;  Budapest,  Copenhagen,  Minsk,  Helsinki,  and  Warsaw,  as 
wind  field  maps  from  the  Taglicher  Wetterbericht  (daily  weather  bulletin)  allowed 
determination  of  wind  velocity  most  readily  for  them.  Resolution  for  these  maps  were 
2.5  km/h  in  wind  speed  and  about  10°  in  direction.  Comparison  of  wind  speeds  will 
follow  the  same  methodology  used  with  the  run  over  the  U.S.  with  separate  paired  t-tests 
for  u  and  v  wind  components. 
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4.  Results  of  OMEGA  Validation 


4. 1  Medium  Range  Forecast  Data  Validation 

Wind  field  predictions  from  September  16-18,  1995,  were  compared  against 
actual  analysis.  Using  the  80  verification  points,  the  direction  was  within  one  compass 
location  (e.g.,  if  the  true  direction  was  west,  this  would  include  west-northwest  and  west- 
southwest)  in  56  of  the  cases.  General  wind  field  patterns  seem  qualitatively  to  predict 
the  winds  accurately.  For  the  complete  list  of  comparisons,  please  refer  to  Appendix  B. 

An  assessment  of  the  magnitudes  was  also  performed,  with  the  components  of  the 
predicted  wind  fields  and  the  analysis  wind  fields  compared  in  a  paired  t-test.  Using  a 
95%  confidence  interval,  the  values  of  u-component  were  not  statistically  different  in  the 
two  sets.  The  mean  value  of  the  analysis  values  was  17.8  and  the  mean  value  of  the 
predicted  values  was  19.0.  The  t-value  was  -0.99.  The  values  of  v-component  were  also 
not  statistically  different.  The  mean  of  the  analysis  values  was  -2.2  and  the  mean  value  of 
the  predicted  values  was  -2.6.  The  t-value  was  0.45.  This  leads  us  to  accept  the  null 
hypothesis  that  the  two  sets  are  the  same.  (Rejection  occurs  with  a  t-value  greater  than 
2.0  or  less  than  -2.0.)  Root-mean-square  errors  for  the  data  were  approximately  10  knots 
for  u-component  and  9  knots  for  v-component.  This  is  good  considering  errors  from  the 
instruments  measuring  these  winds  are  on  the  order  of  6.5  knots. 

Upon  breakdown  of  the  weather  data,  more  bad  data  points  come  from  inclusion 
of  the  station  at  Pittsburgh.  A  separate  paired  t-test  run  without  Pittsburgh  data  showed 
good  comparability.  The  54  data  points  were  compared  using  the  paired  t-test,  resulting 
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in  a  t-value  of  -0.52  and  p-value  of  0.61  for  u-component,  showing  excellent  similarity. 
An  identical  number  of  v-component  data  points  were  compared,  resulting  in  a  t-value  of 
-0.33  and  p-value  of  0.74.  Therefore,  the  indication  is  choosing  stations  along  the 
boundary  cells  produce  worse  results.  This  is  to  be  expected  with  noise  from  the 
boundary  conditions  affecting  the  wind  fields. 

Another  arrangement  of  the  data  was  made  leaving  in  all  three  weather  stations, 
but  removing  the  300  mb  data.  The  60  data  points  were  compared  using  the  paired  t-test, 
resulting  in  a  t-value  of  -0.80  and  p-value  of  0.43  for  u-component.  The  values  for  the  v- 
component  were  -0. 15  and  0.88,  respectively.  We  also  expect  the  higher  wind  field 
assessments  to  contribute  to  error  because  of  larger  velocities  and  proximity  to  the  upper 
boundary  level  affecting  the  domain. 

Figures  4-6  show  the  errors  for  the  wind  field  predictions  based  on  the  output 
from  OMEGA  and  National  Weather  Service  maps.  National  Weather  Service  maps  at 
850  mb  and  OMEGA  generated  wind  fields  are  shown  in  Appendix  C.  The  data  points 
represent  the  National  Weather  Service  observations  with  components  derived  from  the 
maps.  The  difference  between  the  observed  values  and  the  predicted  values  is  plotted. 
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Figure  4.  Error  plot  vs  case  number  of  (a)  u-component  and  (b)  v-component  wind 
speeds  for  Dayton,  Ohio  from  16-18  September,  1995,  beginning  OOZ  at  12-hour 
intervals.  Points  1-7  are  at  300  mb;  8-14  are  at  500  mb;  15-21  are  at  700  mb;  22-28  are  at 
850  mb. 
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4.2  High  Resolution  Analysis  System  Data  Validation 


Wind  field  predictions  from  April  25-28,  1986,  were  compared  against  analysis 
maps  from  the  Taslicher  Wetterbericht.  There  were  60  cases  of  paired  data.  In  only  25 
of  the  cases,  the  predicted  wind  direction  was  within  one  compass  location  of  the 
analysis.  General  wind  field  patterns  seem  to  predict  winds  accurately,  though,  as  much 
complicated  wind  flow  affected  the  area  during  this  time.  For  a  complete  list  of 

comparisons,  please  refer  to  Appendix  B. 

An  assessment  of  magnitudes  was  performed  using  components  of  predicted  wind 
fields  and  analysis  wind  fields  compared  in  a  paired  t-test.  Using  a  95%  confidence 
interval,  the  values  of  the  u-component  were  not  statistically  different  in  the  two  sets  with 
a  t- value  of  -0.79.  The  values  of  the  v-component  were  also  within  the  95%  confidence 
interval,  but  with  a  value  of  1.95,  just  below  the  rejection  value  of  2.0.  The  root-mean- 
square  values  are  also  worse  than  the  case  with  MRF  data,  with  error  values  of  15  knots 
for  u-component  and  1 1  knots  for  v-component. 

Since  the  data  from  28  April  1986  was  suspected  of  being  faulty,  a  separate  test 
was  performed  for  25-27  April  1986,  leaving  45  data  points.  The  u-component  t-value 
was  0.37  with  a  p-value  of  0.71,  excellent  comparison.  The  v-component  t-value  was 
1 .54  with  a  p-value  of  0. 13.  The  rms  error  was  approximately  10  knots  for  u-component 
and  8  knots  for  v-component,  comparable  to  the  results  with  MRF  data. 

On  the  following  page.  Figure  7  shows  the  results  from  the  run  over  Europe.  The 
data  points  represent  the  difference  between  components  derived  from  daily  weather 
maps  and  predicted  values. 
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With  such  drastic  changes  in  wind  flow  in  this  region,  the  resolution  of  the 
domain  itself  would  lend  difficulty  to  predicting  winds  at  a  particular  location  and  time. 
Running  the  model  for  more  than  one  day  continuously  generated  an  increasing  error 
wave  in  the  w-component  winds  in  the  upper  portion  of  the  domain,  which  generally  add 
to  the  error  in  the  domain.  The  predictions  matching  the  wind  direction  were  generally  in 
the  first  two  days  of  the  run. 

The  most  important  factor,  however,  seemed  to  be  the  data  for  28  April  1986. 
Unrealistic  pressures  under  900  mb  at  the  surface  were  predicted  throughout  the  domain, 
and  vertical  wind  speeds  of  over  120  knots  were  predicted.  This  was  the  case  regardless 
of  when  OMEGA  started,  and  casts  suspicion  on  the  data  itself.  These  errors  begin  on 
28  April,  1986,  and  the  model  is  unable  to  recover.  These  errors  were  present  whether  or 
not  a  source  term  was  incorporated  in  the  model. 

Prediction  of  general  wind  flow  patterns  qualitatively  seems  reasonable. 
Deposition  predicted  from  the  Eulerian  tracer  capability  was  compared  to  actual 
deposition  maps  shown  in  Figures  8  and  9,  where  the  maps  for  predicted  deposition  were 
derived  from  the  model.  Appendix  D  shows  the  deposition  maps  actually  produced  by 
OMEGA.  As  mentioned  earlier,  this  provides  a  proof-of-concept  that  OMEGA  could 
reliably  predict  the  deposition,  but  concentrations  are,  unfortunately,  not  available.  An 
assumption  with  the  qualitative  assessment  is  the  dispersion  of  the  tracer  would  generally 
follow  the  same  distribution  but  differ  by  an  unknown  factor. 
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Figure  8.  Actual  deposition  maps  (Smith,  1989:11)  (left)  compared  to  OiVEEGA 
derived  footprint  for  a  1500  meter  release  height 
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4.3  Comment  on  Results 


The  size  of  the  domain  selected  allows  only  a  proof-of-concept  and  not  a  full 
comparison  with  the  ATMES  exercise  participants.  The  smaller  area  selected  makes  the 
boundary  conditions  more  critical  and  provides  OMEGA  with  a  more  stringent  test  of  its 
ability  to  predict  wind  fields. 

Comparison  of  the  wind  fields  using  a  paired  t-test  does  not  reflect  the  standard 
meteorology  test  using  root-mean-square  errors,  but  it  does  indicate  how  OMEGA 
performs  overall  in  prediction  of  the  wind  fields.  Deposition  predictions  depend  less 
upon  the  short-term  variations  in  the  wind  field  than  longterm  projections,  as  seen  by  the 
pattern  for  deposition,  which  compares  well  for  the  first  48  hours  after  the  release.  As  a 
means  to  compare  with  meteorological  models,  values  for  rms  errors  seem  very 
reasonable  for  MRF  data  but  seem  high  for  HIRAS  data 

The  runs  shown  indicate  the  potential  of  OMEGA  to  fulfill  the  needs  of  the 
Defense  Nuclear  Agency.  A  fuller  analysis  at  this  time  is  not  practicable  as  the  model  is 
currently  under  revision.  Some  enhancement  to  the  model  such  as  inclusion  of  a  varying 
source  with  Eulerian  tracers,  capability  to  resume  the  model,  and  capability  to  use  larger 
domains  would  be  required  in  order  to  use  the  same  validation  method  used  in  ATMES. 
As  mentioned  in  the  literature  review,  two  runs  cannot  validate  the  model.  The  model 
would  require  runs  under  various  scenarios  to  permit  the  confidence  in  it  required  for 
assessment  of  limitations  and  accuracy.  These  runs  do  show  the  promise  and  unique 
features  of  OMEGA  at  work. 
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5.  Conclusion 


5. 1  Achievement  of  Objectives 

The  primary  goal  of  the  research  was  to  determine  how  OMEGA  performs 
ingesting  HIRAS  data  and  predicting  deposition  from  a  variable  release  source  term  such 
as  found  in  the  Chernobyl  Nuclear  Accident.  The  intermediate  goal  was  to  describe 
OMEGA’S  performance  with  sample  runs  over  the  United  States  using  MRF  data  by 
comparing  wind  fields  with  weather  data.  These  goals  would  partially  assess  the  validity 
of  OMEGA  to  meet  the  needs  of  the  Defense  Nuclear  Agency. 

The  results  are  favorable.  OMEGA  has  the  capability  to  ingest  data  in  various 
formats  through  its  preprocessor  to  predict  wind  fields  and  meteorological  conditions. 

The  domain  size  and  length  of  simulation  permitted  are  also  adequate  for  a  proof-of- 
concept  validation  of  its  capabilities.  Limitations  with  the  means  of  tracking  particles  do 
not  easily  lend  OMEGA  to  predicting  deposition  with  a  variable  source  term,  but  the 
capability  does  exist  and  use  of  Eulerian  tracer  mode  can  be  used  to  simulate  a  worst-case 
scenario  likely  of  interest  to  the  Defense  Nuclear  Agency.  Based  on  limited  comparison, 
prediction  of  wind  fields  was  within  the  95%  confidence  interval  for  both  domains  and 
data  sets  used. 

5.2  Recommendations 

OMEGA  appears  to  satisfy  the  physical  requirements  for  accurately  predicting 
meteorological  conditions.  More  extensive  study  could  be  taken  in  areas  with  greater 
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numbers  of  weather  stations  in  the  domain  for  greater  statistical  support  for  a  quantitative 
comparison  of  generated  wind  fields  with  analysis  data.  Another  option  would  be  to 
compare  gridded  model  output  with  gridded  analysis  data  to  increase  the  data. 

The  true  need,  however,  is  to  predict  air  concentrations  and  deposition  values.  A 
detailed  look  at  the  Lagrangian  particles  with  their  deviation  values  would  allow 
prediction  of  concentration  and  deposition  by  integrating  over  the  dispersive  range  of 
each  particle  along  each  time  step  (quite  time-consuming).  Once  this  is  performed,  it 
would  assess  the  ability  to  forecast  risk  and  hazard  levels  for  variable  releases,  and  could 
then  be  incorporated  into  the  model’s  framework. 

Improved  user  features  in  the  model  would  also  enhance  future  validation  efforts. 
The  model  is  limited  by  inability  to  resume  after  a  run  has  completed.  Future  versions  of 
the  model  should  allow  restart  capability  by  saving  wind  fields,  particle  locations,  and 
meteorological  conditions.  This  would  allow  greater  flexibility  in  releases  and  allow  the 
most  recent  weather  data  for  the  model.  The  model  could  then  be  allowed  to  run  for  the 
entire  release  period  and  beyond,  permitting  a  true  means  of  comparison  with  models 
participating  in  the  ATMES  exercise. 

Vertical  accelerations  are  also  a  concern  in  the  wind  fields.  The  magnitude  of 
these  accelerations  surpasses  reality  during  extended  runs  (and  sometimes  becomes 
suspect  depending  on  the  given  data.  This  raises  questions  about  the  vertical  layers  and 
how  the  boundary  conditions  at  the  top  of  the  domain  are  resolved.  Further  stabilization 
of  the  boundary  layer  at  the  top  is  necessary  to  prevent  these  unrealistic  vertical 
accelerations. 
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Finally,  it  must  be  remembered  that  OMEGA  is  still  under  revision.  The  research 
was  conducted  in  a  manner  to  highlight  the  capabilities  and  strengths  of  the  model  while 
identifying  concerns  the  Defense  Nuclear  Agency  may  have  regarding  the  model’s 
development.  We  have  seen  great  improvements  and  impressive  capabilities  of  the 
model. 


Appendix  A:  OMEGA  Generated  Grids  for  MRF  and  HIRAS  Validation 


Grid  for  MRF  Validation. 


Grid  for  HIRAS  Validation. 


Appendix  B:  Wind  Field  Tables  for  MRF  and  HIRAS  Validation 


Note;  Case  numbers  are  not  in  the  same  order  as  the  figures  show. 


MIDWEST 


CASE  CITY  DATE  TIME 

MB 

U 

UPRED  V 

VPRED 

1  DAYTON 

9/16/95  OOZ 

300 

14 

42 

6 

-5 

2  DAYTON 

9/16/95  OOZ 

500 

11 

20 

11 

5 

3  DAYTON 

9/16/95  OOZ 

700 

7 

10' 

7 

10 

4  DAYTON 

9/16/95  OOZ 

850 

0 

1 

15 

6 

“  5  DAYTON 

9/16/95  12Z 

300 

80 

51 

0 

-13 

^  6  DAYTON 

9/16/95  12Z 

500 

14 

25 

6 

2 

^  7  DAYTON 

9/16/95  12Z 

700 

14 

11 

6 

7 

8  DAYTON 

9/16/95  12Z 

850 

10 

3 

0 

2 

9  DAYTON 

9/17/95  OOZ 

300 

65 

51 

b 

-13 

10  DAYTON 

9/17/95  OOZ 

500 

14 

25 

-6i 

1 

11  DAYTON 

9/17/95  OOZ 

700 

9 

12 

4 

5 

12  DAYTON 

9/17/95  OOZ 

850 

-6i 

3 

14' 

5 

13  DAYTON 

9/17/95  12Z 

300 

40! 

46! 

Oi 

-4 

14  DAYTON 

9/17/95  12Z 

500! 

28! 

27' 

-11 

-3 

15  DAYTON 

9/17/95  12Z 

700! 

14 

8i 

-14i 

0 

16  DAYTON 

9/17/95  12Z 

850! 

9' 

7' 

-4} 

4 

17  DAYTON 

9/18/95  OOZ 

300! 

23^ 

45! 

-10! 

3 

18  DAYTON 

9/18/95  OOZ 

500! 

14: 

20 

-14i 

-14 

19  DAYTON 

9/18/95  OOZ 

700! 

9! 

-T 

-4! 

-14 

20  DAYTON 

9/18/95  OOZ 

850! 

7' 

-4! 

-71 

4 

21  DAYTON 

9/18/95  12Z 

300! 

32' 

35 1 

-13i 

-13 

22  DAYTON 

9/18/95  12Z 

500! 

23! 

22 

-lOi 

-16 

23  DAYTON 

9/18/95  12Z 

7001 

7' 

6 

-71 

-9 

24  DAYTON 

9/18/95  12Z 

850! 

-4; 

2: 

-9i 

-6 

25  DAYTON 

9/19/95  OOZ 

300! 

28 1 

36 

-28! 

-20 

26  DAYTON 

9/19/95  OOZ 

500! 

14! 

18! 

-141 

-16 

27  DAYTON 

9/19/95  OOZ 

7001 

0! 

5! 

-35! 

-15 

28  DAYTON 

9/19/95  OOZ 

850! 

-9! 

-7! 

-4i 

-7 

29  HUNTINGTON! 

9/16/95  OOZ 

3001 

50 

47' 

01 

-6 

30  HUNTINGTON! 

9/16/95  OOZ 

:  500! 

6! 

191 

14! 

10 

31  HUNTINGTON! 

9/16/95  OOZ 

7001 

11: 

61 

11! 

8 

32  HUNTINGTON! 

9/16/95  OOZ 

850! 

4! 

-1: 

91 

10 

33  HUNTINGTON! 

9/16/95  12Z 

300! 

691 

47' 

-29! 

-17 

34  HUNTINGTON! 

9/16/95  12Z 

500! 

28! 

25  i 

1li 

2 

35  HUNTINGTON! 

9/16/95  12Z 

7001 

4! 

8! 

91 

7 

36  HUNTINGTON! 

9/16/95  12Z 

8501 

71 

3! 

7.11 

8 

37  HUNTINGTON! 

9/17/95  OOZ 

300! 

75  i 

46! 

01 

-1 8 

38!  HUNTINGTON! 

9/17/95  OOZ 

500! 

28! 

241 

12! 

b 

39  HUNTINGTON! 

9/17/95  OOZ 

700! 

18! 

8! 

81 

5 

40  HUNTINGTON! 

9/17/95  OOZ 

850! 

8i 

21 

181 

11 

41  HUNTINGTON! 

9/17/95  12Z 

300! 

401 

541 

Oi 

-4 

42' HUNTINGTON! 

9/17/95  12Z 

500! 

23! 

30! 

-101 

-2 

43;  HUNTINGTON! 

9/17/95  12Z 

700! 

14! 

19! 

-141 

3 

44  HUNTINGTON! 

9/17/95  12Z 

8501 

14! 

lOi 

-141 

6 

45  HUNTINGTON! 

9/18/95  OOZ 

3001 

281 

50! 

-111 

1 

46  HUNTINGTON! 

9/18/95  OOZ 

500! 

14! 

22! 

-14! 

-8 

47  HUNTINGTON! 

9/18/95  OOZ 

700! 

141 

91 

-14! 

-7 

48  HUNTINGTON! 

9/18/95  OOZ 

850! 

9 

2’ 

-4i 

-7 

49  HUNTINGTON! 

9/19/95  OOZ 

500: 

10 

14! 

-23! 

-15 

50  HUNTINGTON! 

9/19/95  OOZ 

700 

0! 

2 

-10! 

-15 

B.l 


MIDWEST 


51  HUNTINGTONI 

9/19/95  OOZ 

850! 

-2: 

•1 

0! 

-7 

52  HUNTINGTON: 

9/19/95  12Z 

300: 

10 

33: 

-23 

-19 

53  PITTSBURGH 

9/16/95  OOZ 

300 

30 

35 

0 

-4 

54  PITTSBURGH 

9/16/95  OOZ 

500. 

20 

24 

0 

0 

55  PITTSBURGH 

9/16/95  OOZ 

700: 

10 

6: 

O' 

2 

56  PITTSBURGH 

9/16/95  OOZ 

850 

8 

2 

19 

5 

57>ITTSBURGH 

9/16/95  12Z 

300; 

30 

40 

0 

-1 

58  PITTSBURGH 

9/16/95  12Z 

500! 

11 

26 

28 

9 

59  PITTSBURGH 

9/16/95  12Z 

700i 

18 

14 

8: 

10 

60  PITTSBURGH 

9/16/95  12Z 

850: 

14! 

3 

14 

5 

eV  PITTSBURGH 

9/17/95  OOZ 

3001 

65 

47 

0! 

-8 

62  PITTSBURGH 

9/17/95  OOZ 

500! 

14! 

26 

6: 

4 

63  PITTSBURGH 

9/17/95  OOZ 

700! 

4i 

15 

9' 

11 

64  PITTSBURGH 

9/17/95  OOZ 

850: 

6 

5: 

14: 

10 

65  PITTSBURGH 

9/17/95  12Z 

300! 

37 

44; 

15! 

-4 

66  PITTSBURGH 

9/17/95  12Z 

1  5001 

201 

231 

01 

-1 

67  PITTSBURGH 

9/17/95  12Z 

700i 

91 

8i 

-4! 

8 

68  PITTSBURGH  ; 

9/17/95  12Z 

:  8501 

01 

6i 

-2! 

19' 

69  PITTSBURGH  , 

9/18/95  OOZ 

:  3001 

23 

401 

-101 

0 

70! PITTSBURGH  : 

9/18/95  OOZ 

5001 

231 

191 

-lOi 

-21 

71  PITTSBURGH 

9/18/95  OOZ 

700! 

141 

12' 

-6i 

-5 

72:  PITTSBURGH: 

9/18/95  OOZ 

8501 

91 

31 

-41 

0 

73:  PITTSBURGH  ^ 

9/18/95  12Z 

300! 

15: 

3T 

Oi 

-5 

74  PITTSBURGH 

9/18/95  12Z 

5001 

15i 

23! 

01 

-9 

75  PITTSBURGH  : 

9/18/95  12Z 

7001 

14! 

9 

-6i 

-4 

76  PITTSBURGH  ; 

9/18/95  12Z 

850! 

9 

91 

-41 

-9 

77  PITTSBURGH  ! 

9/19/95  OOZ 

3001 

32! 

18! 

-131 

-12 

78: PITTSBURGH  : 

9/19/95  OOZ 

500! 

21! 

241 

-211 

-13 

79  PITTSBURGH  ; 

9/19/95  OOZ 

7001 

7! 

14! 

-7! 

-20 

80!  PITTSBURGH  . 

9/19/95  OOZ 

!  8501 

0I 

2! 

-101 

-8 
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ICASE  CITY 

DATE 

MB 

U  UPREC 

V 

VPRED 

Zl 

1  BUDAPEST 

4/25/86 

500 

11 

13 

11 

6 

2  BUDAPEST 

4/25/86 

850; 

19 

17 

19 

6 

'3  BUDAPEST 

4/25/86 

1021 

0 

5 

-5 

1 

"  "4' BUDAPEST 

4/26/86 

500 

38 

21 

0 

0 

5  BUDAPEST 

4/26/86 

850 

22 

10 

0 

0 

6  BUDAPEST 

4/26/86 

1033 

2 

1 

-5 

-2 

7  BUDAPEST 

4/27/86 

500 

16 

12 

16 

9 

8  BUDAPEST 

4/27/86 

850 

11 

20; 

11 

6 

9  BUDAPEST 

4/27/86 

1021 

4. 

11 

4 

2 

10  BUDAPEST 

4/28/86 

500  [ 

8 

-6' 

8 

9 

11  BUDAPEST 

4/28/86 

850^ 

-8 

22 

8 

-7 

12  BUDAPEST 

4/28/86 

896  i 

0 

12: 

-5 

-4 

13  COPENHAGEN 

4/25/86 

500  [ 

O' 

IT 

43 

36 

14  COPENHAGEN 

4/25/86 

8501 

8 

-2: 

8 

8 

15  COPENHAGEN 

4/25/86 

1030! 

-4 

-9i 

4 

0 

16  COPENHAGEN 

4/26/86 

500 

-8 

7: 

20 

23 

17  COPENHAGEN 

4/26/86 

850! 

4 

2: 

-4 

3 

18  COPENHAGEN 

4/26/86 

1037' 

2 

-2! 

-5 

3 

19  COPENHAGEN 

4/27/86 

500' 

10' 

3! 

25 

21 

20  COPENHAGEN 

4/27/86 

8501 

-10! 

-5i 

4 

1 

21  COPENHAGEN 

4/27/86 

10301 

Oi 

-5i 

0; 

-5 

22  COPENHAGEN 

4/28/86 

500! 

-23! 

2! 

23 

28 

23  COPENHAGEN 

4/28/86 

850! 

11' 

17! 

11 

2 

24  COPENHAGEN 

4/28/86 

884  i 

01 

18i 

5' 

1 

25  HELSINKI 

4/25/86 

500! 

23! 

27' 

23 

21 

26  HELSINKI 

:  4/25/86 

850! 

4, 

21 

10, 

5 

27  HELSINKI 

4/25/86 

10491 

-4; 

-7' 

-4 

-2 

281  HELSINKI 

4/26/86 

500! 

-8! 

17! 

20 1 

17 

29' HELSINKI 

4/26/86 

8501 

-19i 

17' 

19' 

2 

30  HELSINKI 

,  4/26/86 

1048! 

-5: 

-1! 

-2: 

10 

31  HELSINKI 

4/27/86 

500  i 

17 

241 

40: 

21 

32  HELSINKI 

4/27/86 

850! 

12! 

10! 

30! 

22 

33  [HELSINKI 

4/27/86 

10351 

01 

3! 

5i 

4 

34i  HELSINKI 

4/28/86 

5001 

161 

27' 

16! 

-9 

35 1  HELSINKI 

4/28/86 

8501 

161 

-421 

16! 

50 

36  HELSINKI 

i  4/28/86 

899! 

Oi 

251 

5i 

-22 

37' MINSK 

;  4/25/86 

500! 

Oi 

-5i 

11: 

10 

38  [MINSK 

i  4/25/86 

850T 

0! 

-2! 

22; 

18 

391  MINSK 

;  4/25/86 

1023! 

-41 

-3! 

4> 

8 

40  MINSK 

4/26/86 

5001 

16! 

-101 

16i 

13 

41  MINSK 

4/26/86 

850! 

-IT 

-71 

11 

21 

42  MINSK 

i  4/26/86 

1032! 

-5i 

-4i 

0! 

9 

43  MINSK 

4/27/86 

5001 

16! 

-7! 

16! 

-3 

44  MINSK 

:  4/27/86 

8501 

Oi 

-18! 

22: 

22 

45 1  MINSK 

4/27/86 

1017! 

-4i 

-7! 

-4i 

9 

46  [MINSK 

4/28/86 

500! 

32: 

27' 

0; 

-9 

47’ MINSK 

4/28/86 

850! 

-16 

-17! 

16i 

10 

48  MINSK 

4/28/86 

891 

-5 

-9' 

-2 

8 
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49  WARSAW 

4/25/86  500 

16 

9: 

16 

11 

50  WARSAW 

4/25/86  850 

0- 

6 

11 

15 

51  WARSAW 

4/25/86  1014 

-4 

0 

-4 

8 

52  WARSAW 

4/26/86  500 

16 

12 

16 

11 

'  "  53  WARSAW 

4/26/86  850 

16 

4 

16 

12 

54  WARSAW 

4/26/86  1032 

-5 

3 

2 

0 

■  “  55  WARSAW 

4/27/86  500; 

16 

11 

16 

0 

56‘ WARSAW 

4/27/86  850 

-8 

"-lO" 

20 

5 

57  WARSAW 

4/27/86  10i5 

-6 

- 

"  -4 . 

4 

58  WARSAW 

4/28/86  500; 

-1l" 

-3 

. 59  WARSAW 

4/28/86  850 

35 

-11 

-28 

60  WARSAW 

4/28/86  895 

-4 

18 

4 

-15 
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Appendix  C:  National  Weather  Service  Maps  and  OMEGA  Generated  Maps  for 


MRF  Validation 


09/ 16  OOZ  850hPa  Analysis  | 

160200Z  Scfvere  Tstm  otlk  valid  until  161200Z  (SHOW  REGION) 


C.2 


OMIilGA  Wind  Field  i«n  09/16  OOZ  850  ml) 


C.3 


C.4 


OMEGA  Wind  Field  Idr  ()9/.l6  I2Z  85(1  mb 


C.5 


C.7 


OMEGA  Wind  Field  for  09/17  12Z  850  mb 


522 


CIO 


OMEGA  Wind  Field  for  09/18  OOZ  850  mb 


09/18 12Z  850hPa  Analysis  | 

181200Z  Severe  Tstm  oUk  valid  until  191200Z  (SHOW  REGION). 
SviTstim  Watch  #WS904  valid  unUlDC^  I 


C.ll 


C.12 


OMEGA  Wind  Field  for  09/18  I2Z  850  mb 


C.13 


OMEGA  Wind  Field  lor  09/19  OOZ  850  mb 


Appendix  D:  Sample  Maps  for  Deposition  Using  HIRAS  Data 
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OMEGA  Generated  Footprint  for  Eulerian  Tracer  at  1500  meters,  04/27/86  OOZ. 


D.3 


OMEGA  Generated  Footprint  for  Eulerian  Tracer  at  1500  meters,  04/27/86  12Z. 


D.4 


OMEGA  Generated  Concentration  at  300  mb  for  Eulerian  Tracer  at  1500  meters, 
04/28/86  OOZ. 


Appendix  E:  Chernobyl  Source  Term  Defined  by  Fortran  90  Code 


c  Subroutine  Chernobyl  to  model  nuclear  accident  on  26  April  1986, 
c  OOOOZ.  Source  term  is  modeled  according  to  the  Technical 
c  Specifications  Document  in  "Evaluation  of  Long  Range  Atmospheric 
c  Transport  Models  using  Environmental  Radioactivity  Data  From  the 
c  Chernobyl  Accident",  the  ATMES  Report, 
c 

program  chemobyl 
parameter(numsrc=  1 2) 
integer  i,  maxid,  iprt_idO(numsrc) 

real  timeinj,  prt_tim(numsrc+l),  prt_mass(numsrc),  prt_dO(numsrc), 
&  prt_lonO(numsrc),  prt_IatO(numsrc),  prt_altO(numsrc), 

&  prt_rhdO(numsrc),  prt_vapO(numsrc),  prt_sxO(numsrc), 

&  prt_syO(numsrc),  prt_szO(numsrc),  fact 
character*6  modeadm 

fact  =  4.5E+07*2.29E-25 
timeinj  =  3600. 
modeadm  =  'tracer' 

prt_tim(l)  =  86400. 
pn_tim(2)  =  86400+6.*60.*60. 
prt_mass(l)  =  2.2E+16*fact*.2/6. 
prt_mass(2)  =  2.2E+16*fact*.8/18. 
prt_mass(3)  =  7.0E+15*fact/24. 
prT_mass(4)  =  5.51E+15*fact/24. 
prt_mass(5)  =  4.1  lE+15*fact/24. 
prt_mass(6)  =  3.01E+15*fact/24. 
prt_mass(7)  =  prt_mass(6) 
prt_mass(8)  =  prt_mass(4) 
prt_mass(9)  =  6.31E+15*fact/24. 
prt_mass(10)  =  8.11E+15*fact/24. 
prt_mass(l  1)  =  8.91E+15*fact/24. 
prt_mass(12)  =  l.llE+14*fact/24. 
prt_altO(l)  =  1500. 

do  i=l,  numsrc 

prt_d0(i)=l.e-20 
iprt_idO(i)=i 
prt_lonO(i)  =  30.15 
prt_latO(i)  =  51.17 
prt_rhdO(i)=  1  .e-20 
prt_vapO(i)=  1  .e-20 
prt_sx0(i)=400. 
prt_sy0(i)=400. 
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prt_sz0(i)=400. 

c  Now  set  altitude,  each  variation. 

if  (i  .gt.  1  .and.  i  .It.  4)  prt_altO(i)  =  600. 
if  (i  .gt.  3)  prt_altO(i)  =  300. 

c  Set  time  of  release,  each  variation.  ' 

if  (i  .gt.  2)  prt_tim(i)  =  (i-l)*86400. 
enddo 

prt_tim(13)  =  prt_tim(12)+86400. 

c  Now  let's  write  the  output  to  a  file 

open(unit=l,  file='case.adm',  status='unknown’) 
write  (1,  *)  maxid 
write  (1,  *)  numsrc 
write  (1,100)  timeinj 
100  format(f6.0) 

write  (1,  110)  modeadm 
110  format(a6) 

do  i=l,  numsrc 

write  (1,  120)  iprt_id0(i),  prt_lon0(i),  prt_lat0(i),  prt_alt0(i), 
&  prt_tim(i),  prt_tim(i+l),  prt_d0(i),  prt_mass(i),  prt_rhd0(i), 

&  prt_vap0(i),  prt_sx0(i),  prt_sy0(i),  prt_sz0(i) 

120  format(i8,12el5.7) 
enddo 
stop 
end 
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